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ABSTRACT 

We present the results of a 5 GHz survey with the Very Large Array (VLA) and the expanded VLA, 
designed to search for short-lived 1 day) transients and to characterize the variability of radio 
sources at milli-Jansky levels. A total sky area of 2.66 deg 2 , spread over 141 fields at low Galactic 
latitudes (b = 6-8 deg) was observed 16 times with a cadence that was chosen to sample timescales 
of days, months and years. Most of the data were reduced, analyzed and searched for transients in 
near real time. Interesting candidates were followed up using visible light telescopes (typical delays of 
1-2 hr) and the X-Ray Telescope on board the Swift satellite. The final processing of the data revealed 
a single possible transient with a flux density of f v = 2.4 mJy. This implies a transients sky surface 
density of > 1.8 mJy) = O.OSQlo Jg^t-choss deg~ 2 (1, 2-cr confidence errors). This areal density 
is consistent with the sky surface density of transients from the Bower et al. survey extrapolated to 
1.8 mJy. Our observed transient areal density is consistent with a Neutron Stars (NSs) origin for these 
events. Furthermore, we use the data to measure the sources variability on days to years time scales, 
and we present the variability structure function of 5 GHz sources. The mean structure function shows 
a fast increase on ps 1 day time scale, followed by a slower increase on time scales of up to 10 days. 
On time scales between 10-60 days the structure function is roughly constant. We find that <; 30% of 
the unresolved sources brighter than 1.8 mJy are variable at the > 4-er confidence level, presumably 
due mainly to refractive scintillation. 

Subject headings: radio continuum: general — stars: neutron — techniques: photometric 



1. INTRODUCTION 

Radio surveys of the sky in the time domain have of- 
ten been used to identify new astrophysical phenomena. 
Highly variable radio sources can serve as signposts to 
compact, high energy objects which are accompanied by 
high magnetic fields and/or relativistic particle acceler- 
ation. Radio variability from quasars and 7-ray bursts 
(Dent 1965; Frail et al. 1997) was used to infer bulk 
relativistic motions in these objects (Rees 1967; Good- 
man et al. 1987). Notable new phenomena identified 
from radio time-domain surveys include the discovery of 
the first pulsars (Hewish et al. 1968), the Galactic high- 
energy binary LSI+61°303 (Gregory and Taylor 1978), 
the anomalous variability of 4C 21.53 that lead to the 
discovery of millisecond pulsars (Backer et al. 1982), and 
the still-mysterious extreme scattering events (Fiedler et 
al. 1987). 

More recent surveys have found several new types of 
radio transients whose identity has remained unknown or 
not well understood (e.g., Hyman et al. 2005; McLaugh- 
lin et al. 2006; Bower et al. 2007; Lorimer et al. 2007; 
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Niinuma et al. 2007; Kida et al. 2008; Matsumura et al. 
2009). 

Specifically, Bower et al. (2007) re-analyzed 944 epochs 
of Very Large Arraj0 (VLA) observations, taken about 
once per week for twenty two years, of a single calibra- 
tion field. These authors discovered a total of ten tran- 
sients, eight in the 5-GHz band and two in the 8-GHz 
band. Eight of these transients were detected in a single 
epoch. Therefore, their duration is shorter than the time 
between successive epochs (one week) and longer than 
the exposure time (20min). Moreover, the majority of 
these sources do not have any optical counterpart coin- 
ciding with their position. The lack of optical counter- 
parts down to limiting magnitudes of 27.6 in (/-band and 
26.5 in i?-band is especially puzzling and significantly 
limits the classes of objects that can be associated with 
these events (Ofek et al. 2010). 

In a possibly related work, Kuniyoshi et al. (2006), 
Niinuma et al. (2007), and Kida et al. (2008) reported 
a search for radio transients using an East- West inter- 
ferometer of the Nasu Pulsar Observatory (located in 
Tochigi Prefecture, Japan) of Waseda University. To 
date, this program reported 11 bright radio transients 
with flux densities above 1 Jy in the 1.4-GHz band. 

Recently, in Ofek et al. (2010) we suggested that the 
properties of the single epoch "Bower et al. transients" 
and the Nasu transients are consistent with emerging 
from a single class of objects, namely isolated old Neu- 
tron Stars (NS). Specifically, the NS hypothesis is consis- 

8 The Very Large Array is operated by the National Radio As- 
tronomy Observatory (NRAO), a facility of the National Science 
Foundation operated under cooperative agreement by Associated 
Universities, Inc. 
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tent with the rate, energetics, sky surface density, source 
number count function and the lack of optical counter- 
parts. 

In this paper we present a new VLA survey for radio 
transients and variables at low Galactic latitudes. Our 
main motivation for this survey was to detect more ex- 
amples of this new class of short-lived radio transients, 
with the goal of identifying them in real-time in order to 
find their counterparts at other wavelengths for further 
study. A second, and equally important motivation for 
this survey was to characterize the transient and variable 
radio sky with a sensitivity and cadence which had not 
been carried out previously. 

The organization of this paper is as follows. In $2] we 
provide a summary of previous radio transient and vari- 
ability surveys. In $3] we present the observations, while 
the data reduction is outlined in $4] The results from 
our real time transients search are provided in Sj5] <J6] 
present the source catalogs generated in the post survey 
phase. The final post survey transient search is described 
in Sj7] while the sources variability study is presented in 
Sj8] The implications of this study are discussed in Sj9]and 
we summarize in £1101 In addition three appendices dis- 
cussing: flux calibration; the statistics of max/min of a 
time series; and transient areal density calculation in the 
case of a beam with non-uniform sensitivity are provided. 

2. PREVIOUS GHZ SURVEYS FOR TRANSIENTS 
AND VARIABLES 

Existing 0.8-8 GHz surveys have already explored, to 
some extent, the dynamic radio sky with a wide range 
of sensitivities, angular resolution and cadences. How- 
ever, compared with synoptic surveys at higher frequen- 
cies (infra-red to 7-rays) the radio sky remains poorly 
explored. In Table Q] we summarize past synoptic ra- 
dio surveys. For each survey we list also the number of 
transients, as well as variables which vary by more than 
50%. We note however, that comparison of these num- 
bers is complicated due to several factors. A radio image 
may be accomplished either through a single pointing, 
or adding several scans taken at different times. If the 
time span, St, containing all the observation composing 
a single "epoch" is larger than the transient duration (or 
variability time scale) then the survey sensitivity to tran- 
sients (variables) is degraded. Additionally, the proba- 
bility of detecting significant variability depends on St 
and the typical time scale between epochs (At), through 
the variability structure function. Depending on the sta- 
tistical method used to define the variability amplitude, 
it may also affected by the number of epochs (N ep ) in 
the survey. Here we provide a summary of some 

of the sky surveys listed in Table Q] 

2.1. Surveys at frequencies below 2 GHz 

Carilli et al. (2003) used a deep, single VLA pointing 
at 1.4 GHz toward the Lockman hole. They found that 
only a small fraction, < 2%, of radio sources above a flux 
density limit of 0.1 mJy are highly (>50%) variable on 19 
day and 17 month timcscales. No transients were iden- 
tified. Frail et al. (1994) imaged a much larger field at 
1.4 GHz toward a 7-ray burst with the Dominion Radio 
Astrophysical Observatory synthesis telescope, making 
daily measurements for two weeks and then on several 



single epochs for up to three months. No transients were 
identified on these timescales, and no sources above a 
flux density limit of 3.5 mJy were seen to vary by more 
than 4<t. 

There have also been a number of wide field surveys 
of the sky at 1.4 GHz. de Vries et al. (2004) used the 
VLA to image a region, toward the South Galactic cap, 
twice on a seven year timescale. No transients were found 
above a limit of 2mJy. Croft et al. (2010; 2011) pre- 
sented results from the Allen Telescope Array Twenty- 
centimeter Survey (ATATS). They surveyed 690 deg 2 of 
an extragalactic field on 12 epochs. They compared 
the individual images with their combined image and 
the combined image with the NRAO VLA Sky Survey 
(NVSS; Condon et al. 1998). No transients were found 
above a flux density limit of 40 mJy in the combined 
image, with respect to the NVSS survey (Croft et al. 
2010). In addition, no transients were found in the indi- 
vidual epochs above flux density of about 100 mJy (Croft 
et al. 2011). 

A systematic search for transients was made between 
the two largest radio sky surveys, Faint Images of the 
Radio Sky at Twenty- Centimeters (FIRST; Becker et al. 
1995) and NVSS (Condon et al. 1998) by Levinson et 
al. (2002). Nine transient candidates were identified. 
Follow-up observations of these established that only one 
was a genuine transient - a likely radio supernova in 
NGC4216 (Gal- Yam et al. 2006; Ofek et al. 2010). We 
note that each FIRST and NVSS image is composed of 
ss 4 overlapping beams of adjacent regions taken typi- 
cally with St ^days (Becker et al. 1995; Condon et al. 
1998; Ofek & Frail 2011). Therefore, if the duration 
(idur) of the Bower et al. transients is shorter than this 
typical time between images composing "one epoch" then 
the sensitivity of these surveys for transients is degraded 
by rj \/4. However, the Levinson et al. survey used a 
flux density limit of 6 mJy which is (<; 2) higher than the 
flux limit of these surveys. Therefore, its efficiency for 
Bower et al. like transients is not degraded. 

Bright (> Uy), short-lived transients, have been re- 
ported by the Nasu 1.4 GHz survey (e.g., Matsumura et 
al. 2009; see JQ). Croft et al. (2010; 2011) argued that 
these transients are not real because their implied event 
rate cannot be reconciled with their own survey unless 
this population has sharp cut off at flux densities below 
1 Jy. We note that Croft et al. adopted the Nasu tran- 
sients areal density reported in Matsumura et al. (2009). 
However, this areal density is inconsistent with the rate 
reported in Kida et al. (2008) which is roughly two orders 
of magnitude lower. Furthermore, based on the Nasu sur- 
vey parameters reported in Matsumura et al. (2009) we 
estimate that the areal density of the Nasu transients is 
roughly two orders of magnitude lower than that stated 
in their papei0. 

Here we present an estimate of the rate of transients 
from the Nasu observations: Matsumura et al. (2009) 
reported that they discovered (at that time) nine tran- 
sients over a period of two years (730 days). Assum- 
ing that they have four pairs of antennas (N ant ) each 

9 This inconsistency was already mentioned by Bower & Saul 
(2011). However, they reach a different conclusion than we do. 
More information about the Nasu survey is needed in order to 
resolve this issue. 
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TABLE 1 

Previous GHz Transient and variability Surveys 



V 


Area 


Direction 




N ep 


St 


Af 


rms 


Sources 


Tran. 


Var. 






Ref. 


GHz 


deg 


deg 








mjy 














0.84 


2776 


8 < -30 


~ 45 


2 a 


12 hr 


1 day-20 yr 


2.8 


29730 


15 


~ 10 


14] 


1.4 


0.22 


I = 150, b = +53 


4.5 


3 


6 hrs 


19d, 17m 


0.015 







2% 


1 






1.1 


2.6 


I = 151 b = +24 


60 


16 


12 hrs 


1-12 d 1-3 m 


0.7 


245 


o 


~ 1% 


2 






1.4 


120 


S. Galactic Cap 


5 


2 




7 yr 


0.15 


9086 


o 


1.4% 


3 






1.4 


2500 


b ^ 30 


45 


2 


days 




0.45 


7181 


I 




5,f 


,7] 


1.4 


2870 b 


4-^9 ^ A ^> -1-49 


24' X 2.4' 


~ 1000 


4 linn 


1 d 


300 








8,9,10,11] 


1.4 


0.2 


/ — 57 h — 4-81 


20 


1852 


minutes 


1 day — 23 yr 


2 


10 


o 




19] 


1.4 


690 


2 = 70, b = +64 


150 


2 


months 


15 yr 


3.94 


4408 





~ 0.1% 


4] 




1.4 


690 


2 = 70, b = +64 


150 


12 


> lday 


days— months 


38 


4408 





£ 0.5% 


20 




1.4 


0.2 


phase calib. 




151 


5 min 


days- years 


~ 1 









21 




3.1 


10 


1 = 57, 6 = +67 


100 


2 


months 


15 yr 


0.25 


425 


l c 




12 




4.9 


0.1 


phase calib. 




~ 390 d 


5 min 


days- years 


~ 1 









21 




4.9 


0.69 


Extragalactic 


0.5-15 


2 


60 min 


l-100d 


0.05 









15 




4.9 


23.2 


\b\ < 0.4 


5 


3 


90s 


2 m— 15 yr 


0.2 


2700 





15 


16 




4.9 


500 


|6|<2 


180 


16 


2 min 


1 day-5 yr 


1.6 


1274 


1 


~ 0.5 


18 




4.9 


19924 


75 > <5 > 


210 


2 


--'week 


lyr 


5 


75162 





> 40 


17 




4.9 


0.07 


I = 115, 6 = +36 


5 


626 


20 mim 


1 week-22 yr 


0.05 


8 


-c 





13 




8.5 


0.02 


/ = 115, b = +36 


3 


599 


20 min 


1 week-22y 


0.05 


4 


l f 





13 




8.5 


0.04 


phase calib. 




~ 308s 


5 min 


days- years 


~ 1 









21 




4.9 


2.6 


|6|«7 


15 


16 


50 s 


1 d-2 yr 


0.15 


~ 200 


1 


0.3 - 30% This paper 



Note. — Columns description: AB is the beam full width at half power; 7V ep is the number of epochs; 5t is the time span over which each epoch was 
obtained (see text); Af is the range of time separations between epochs; Sources is the number of persistent sources detected; Tran. is the number of 
transients found by the survey; Var. is the number or percentage of variables showing variability larger than 50%. We note that strong variables arc defined 
differently by each survey. Therefore, these numbers provide only a qualitative comparison between the surveys. References: [1] Carilli ct al. (2003), [2] Frail 
ct al. (1994), [3] dc Vrics ct al. (2004), [4] Croft ct al. (2010), [5] Lcvinson ct al. (2002), [6] Gal- Yam ct al. (2006), [7] Ofck ct al. (2010), [8] Matsumura ct 
al. (2009), [9] Kida ct al. (2008), [10] Kuniyoshi et al. (2007), [11] Matsumura et al. (2007), [12] Bower ct al. (2010), [13] Bower et al. (2007), [14] Bannister 
ct al. (2010), [15] Frail et al. (2003), [16] Becker et al. (2010), [17] Scott (1996), [18] Gregory & Taylor (1986), [19] Bower & Saul (2011), [20] Croft et al. 
(2011), [21] Bell ct al. (2011). 

a Smaller fraction of the sky was observed more than twice. 

The total surveyed area is about 2870 deg 2 , but about 460 deg 2 was surveyed every day. These parameters are deduced from the Kida ct al. (2008) and 
Matsumura ct al. (2009) papers (sec j2j. 
c Marginal detection (4.3cr). Ignored in Figure [Tl 

Mean number of epochs per field - Seven fields were observed on 2732 epochs. 
° In addition, one transient was found by combining two month worth of data and no transients were found by combining lyr worth of data. 
^ In addition, one transient was found by combining two month worth of data and no transients were found by combining lyr worth of data. 
g Mean number of epochs per field - Seven fields were observed on 2154 epochs. 



looking at a different position (near the local zenith 
S w 37 deg) with a field of view of W = 0.4 deg and 
scanning the sky at the sidereal rate, the total sky area 
scanned by their system after two years is « 3.4 x 10 5 deg 2 
(= N ant W x 360cos(37°) x 730). Therefore, the total 
areal density of their transients is the number of tran- 
sients divided by the total scanned sky area and divided 
by 1.61, which is « 1.7 x 10~ 5 deg" 2 . Where the 1.61 
factor is due to the fact that their sensitivity is not uni- 
form within the beam and, and is calculated using Equa- 
tion IC4I in Appendix [D] Assuming the transients du- 
ration is lday than their rate is « 0.02 deg" 2 yr" 1 . We 
note that our derived Nasu transients rate is roughly con- 
sistent with the upper limit on the rate given by Kida et 
al. (2008). Therefore, we speculate that there was a 
confusion between areal density and transient rate in the 
Matsumura et al. (2009) paper. We conclude that if our 
estimate for the areal density for the Nasu transients is 
correct than the ATATS survey results cannot decisively 
rule out the reality of the Nasu transients. 

A comprehensive survey at 0.8 GHz was reported by 
Bannister et al. (2010). They surveyed 2776 deg 2 south 
of S = —30° over a 22-year period. Out of about 
30,000 sources they identified 53 variables and 15 tran- 
sient sources. Recently, Bower & Saul (2011) reported 
a transient search in the fields of the VLA calibrators 
in which no transients were found (Summarized in Ta- 



ble [T] and Figure [lj . Another related work by Bell et al. 
(2011) searched for radio transients in the fields of the 
VLA phase calibrators at 1.4 GHz, 4.9 GHz and 8.5 GHz. 
Based on their survey parameters (Table Q} we estimate 
that their 95% confidence surface density upper limit 
on transients brighter than 8mJy in 1.4 GHz, 4.9 GHz 
and 8.5GHz are 0.19deg~ 2 , 0.13deg~ 2 and 0.50deg" 2 , 
respectively. These values are corrected for the beam 
non-uniformity factor (of 1.61), mentioned earlier. We 
assumed that Bell et al. (2011) searched for transients 
within the full width at half power of the beam. For 
clarity purposes, in Figure [T] we show only the 4.9 GHz 
limit. 

2.2. Surveys at frequencies above 2 GHz 

There have been several additional transient surveys 
carried out at frequencies above 1.4 GHz. At 3.1 GHz 
Bower et al. (2010) report a marginal detection of one 
possible transient (4.3<r) in a 10 deg 2 survey of the Bootes 
extragalactic field. In a five year catalog of radio after- 
glow observations of 75 7-ray bursts, Frail et al. (2003) 
found several strong variables at 5 and 8.5 GHz, but 
no new transients apart from the radio afterglows them- 
selves. 

Two surveys at 5 GHz have specifically targeted the 
Galactic plane. Taylor and Gregory (1983) and Gregory 
and Taylor (1986) used the NRAO 91-m telescope to im- 
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age an approximately 500 deg 2 region from Galactic lon- 
gitude I = 40 deg to / = 220 deg with Galactic latitude 
|6| < 2 deg in 16 epochs over a 5-year period. They identi- 
fied one transient candidate which underwent a Uy flare 
but for which follow-up VLA observations showed no qui- 
escent radio counterpart (Tsutsumi et al. 1995). They 
also claimed tentative evidence for a separate Galactic 
population of strong variables comprising 2% of their 
sources. Support for this comes from Galactic survey of 
Becker et al. (2010) , who find about one half of their vari- 
able source sample (17/39), or 3% of all radio sources in 
the Galactic plane, undergo strong variability on 1-year 
and 15-year baselines. We note that the surface den- 
sity of radio sources in the Galactic plane is only slightly 
higher (w 20%) than at high Galactic latitudes (Helfand 
et al. 2006; Murphy et al. 2007). 

One of the largest variability survey of its kind was 
carried out using the 7-beam receiver on the NRAO 91- 
m telescope (Scott 1996; Gregory et al. 2001). The 
sky from 0° < <5 < 75° was surveyed over two 1-month 
periods in 1986 November and 1987 October (Condon, 
Broderick & Seielstad 1989; Becker et al. 1991; Condon 
et al. 1994; Gregory et al. 1996). The final catalog, 
made by combining both the 1996 and 1997 epochs, con- 
tained 75,162 discrete sources with flux densities > 18 
mjy. Long term variability information was available for 
the majority of the sources by comparing the mean flux 
densities between the 1986 and 1987 epochs. Scott (1996) 
carried out a preliminary analysis of the long-term mea- 
surements and identified 146 highly variable sources, or 
<1% of the cataloged radio sources. 

Eight possible transients in the Scott (1996; Table 5.1) 
list appear in either 1986 or 1987 but are undetected 
in the other epoch (< 2a). Two sources are previ- 
ously identified variables from the Gregory and Taylor 
(1986) survey, while six are flagged as possible false pos- 
itives due to confusion by nearby bright sources. One 
source (B150958.3+103541) was 9 ± 6mJy in 1986 and 
75±7mJy in 1987 but it is in both the FIRST and NVSS 
source catalogs. There are therefore no long-term tran- 
sients identified in the Scott (1996) survey. 

In order to estimate the flux limit above which the 
Scott (1996) comparison between the 1986 and 1987 sur- 
veys is complete, we compared the source numbers near 
the celestial equator, as a function of flux in the two 
publicly available catalogs from 1987 and the combined 
1986/1987 catalog. We found that at flux densities be- 
low about 40mJy the number of sources, as a function 
of flux, in the 1987 catalog is rising slower than that for 
the one of the deeper combined catalog. Therefore, we 
estimate that near the terrestrial equator, the 1987 cat- 
alog is complete above a flux density of about 40 mjy. 
However, these catalogs were made from observations in 
which each point on the sky was observed « 4/cos(<5) 
times taken within a few days. This degrades the sen- 
sitivity of the comparison carried out by Scott (1996), 
for <> 1 day transients, by about \/\. Therefore, we con- 
clude that the Scott (1996) survey is sensitive to short 
term (<S 1 day) transients brighter than about 80mJy 
(= 40v^)- Finally, assuming Scott (1996) did not find 
any transients in two epochs, we put a 2-er upper limit 
on the areal density of <J 1 day transients brighter than 
80mJy, of 9.5 x 10~ 5 deg~ 2 . 



TABLE 2 

List of survey pointings 



Field Name 


RA 


Dec 


A cp 






deg 


deg 


1851- 


-1327 


282.94699 


-13.45500 


16 


1852- 


-1309 


283.09963 


-13.15668 


16 


1853- 


-1233 


283.40439 


-12.56005 


16 


1853- 


-1251 


283.25209 


-12.85836 


16 


1853- 


-1318 


283.40625 


-13.30508 


16 



Note. — List of all 141 fields that were ob- 
served as part of this survey. The number of 
epochs per pointing is marked in N cp . This table 
is published in its entirety in the electronic edi- 
tion of the Astrophysical Journal. A portion of 
the full table is shown here for guidance regarding 
its form and content. 



In summary, despite the heterogeneous nature of these 
GHz surveys, it is clear that the radio sky is relatively 
quiet compared with the 7-rays sky. The fraction of 
strong variables among the persistent radio source pop- 
ulation is 0.1-3% from flux densities of 0.1 mjy to Uy. 
However, the exact percentage of strong variables is still 
uncertain because of the different criteria used by vari- 
ous surveys. We note that within this flux density range, 
radio source populations are dominated by AGN roughly 
above 1 mjy and star forming galaxies dominates the 
source counts below 1 mjy (Condon 1984; Windhorst 
1985). 

The transient areal densities detected by these various 
surveys, as well as our survey, are shown graphically in 
Figure [TJ Also shown in this figure are the persistent 
sources areal densities at different frequencies. This plot 
is further discussed in Sj§] 

3. SURVEY OBSERVATIONS 

We designed a survey to look for transients and vari- 
able sources near the Galactic plane, with typical time 
scales of days to two years at milli-Jansky flux levels. 
We were specially interested in finding "Bower et al. 
transients" , conduct multi- wavelength follow up of these 
events and finding counterparts and studying their spec- 
tral evolution. 

3.1. Survey Design 

We used the VLA to observe 141 pointings along the 
Galactic plane. In order to minimize telescope mo- 
tions we selected all the pointings in four regions. The 
median longitude (I) and latitude (b) of the four re- 
gions are: I = 22.6 deg, b = —6.7 deg; I = 56.6 deg, 
b = -5.5deg; I = 89.7deg, b = -7.8deg; I = 106.0deg, 
b = —6.5 deg. Within each region we selected 26- 
42 pointings within 2.3 deg from the median position 
of each region. Each pointing was selected to have 
no NVSS sources brighter than 1 Jy within 3 deg, no 
NVSS sources brighter than 300 mjy within 1 deg, and no 
source brighter than 100 mjy within the field of view as 
defined by the half power radius. We also rejected fields 
which distance from known Galactic supernova remnants 
(SNR; Green 2001) is within twice the diameter of the 
SNR. The typical distance between pointings in each re- 
gion is about 20'. The final 141 pointings are listed in 
Table H 
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Fig. 1. — Cumulative areal density of radio sources and transients as a function of flux density for various surveys. Different colors 
represent different frequencies as specified in the legend. 95% confidence upper limits from various transient surveys are shown as right- 
angle corners, while measured areal densities are marked as filled circles. All the error bars represent 2-<r confidence intervals. We note that 
the Nasu survey rate is based on our estimate using the survey description in their papers (see {J2}- Because some of the parameters of this 
survey are unknown to us, we increase the error bars for this survey to include a factor of two uncertainty. For the Levinson et al. survey 
we mark the areal density of the single transient found in this survey (a supernova in NGC 42f 6) and also the 95% confidence upper limit 
assuming there are no Bower et al. transients in the survey. For the de Vries et al. (2004) search we use a flux limit of 2 mJy since they 
used the FIRST survey in which each epoch is composed of about four observations of the same field taken within St ssdays. Therefore, 
this degrade their survey flux limit sensitivity to Bower et al. transients by a factor of about \/\. The right-hand-side y axis shows the 
transients rate assuming a transient duration of 0.5 days. Some surveys are excluded from this plot. For example, Becker et al. (20f0) 
restricted their catalog to sources detected in at least two out of three epochs, or that have a confirmed detection at 1.4 GHz. Therefore, 
such surveys are not included here. Also shown are the areal densities of persistent sources at 1.4 GHz and 4.9 GHz based on the FIRST 
and GB87 surveys, respectively (dashed lines). 



3.2. Observations 

These 141 fields were observed on 11 epochs using the 
VLA in the Summer of 2008 and on five epochs using 
the Expanded VLA (EVLA) during the Summer of 2010. 
All observations were made in the compact D configu- 
ration. For the 2008 observations, we added together 
two adjacent 50 MHz bandwidths centered at 4835 and 
4885 MHz with full polarization. For the 2010 observa- 
tions we added together two adjacent 128 MHz sub bands 
centered at 4896 and 5024 MHz with full polarization. 

In 2008 care was taken to ensure that the local sidereal 
start time was the same for each 3-hr epoch (20:30 LST). 
Therefore, each field was observed at the same hour an- 
gle and subsequently the synthesized beam stayed the 
same for each epoch, varying only when antennas were 
taken out of the array. The 2010 observations were taken 
during EVLA shared-risk science commissioning, and so 
some scans lost due to correlator errors and the last two 
epochs began one hour earlier than our 2008 local side- 
real start time. 



We integrated each pointing for about 50 s on average. 
The maximum integration time was 58.5 s and the mini- 
mum was 43.3 s. Additionally, during each 3-hr observing 
run we carried out all necessary calibrations. Amplitude 
calibration was achieved with observations of 3C 286 and 
3C 147 at the start and end of each epoch, respectively. 
Phase calibration was checked every 20-25 min by switch- 
ing to a bright point source within a few degrees of the 
targeted region. We used the following four phase cal- 
ibrators (one per each region): J1911— 201, J1925+211, 
J2202+422, and J2343+538. The total calibration and 
antenna move-time overhead was about 30% of the ob- 
serving time. This overhead on move time could have 
been lowered had we used the fast slew methods from the 
NVSS and FIRST surveys (Condon et al. 1998, Becker 
et al. 1995) with a resulting increase in the number of 
square degrees of sky surveyed per hour. However, since 
we recently found that this method could introduce spu- 
rious transients (Ofek et al. 2010), we adopted a less 
efficient but more robust observing method. 
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In Table [3] we list the time of the UTC midpoint 
for each epoch with some additional information. The 
shortest variability timescale sampled was 24 hr and the 
longest was = 2yr. By design, the cadence of the 2008 
survey was chosen to probe variability timescales be- 
tween a day and a month. A longer 2 yr timescale was 
also sampled by comparing deep images made from the 
2008 and 2010 campaigns (see below). 

4. DATA REDUCTION AND CALIBRATION 

In 2008 the uv data was streamed directly to a disk in 
real time, and a pipeline was ran after each one of the 
four regions were observed. We used the data reduction 
pipeline provided in the Astronomical Image Processing 
System (AIPS) packagcFl For each epoch, the pipeline 
first flagged and calibrated the uv data. It then imaged a 
30-arcmin-wide field around all 141 pointings, deconvolv- 
ing down to three times the rms noise, and restoring the 
image with a robust weighted beam. No self-calibration 
was done. The VLA data rates (« 30 Mbytes hr -1 ) 
and the D-configuration image requirements (512 pix- 
els, 3.6" pixel -1 ) were so modest that the entire pipeline 
reduction and the variable source analysis (see was 
completed before the VLA finished observing the next 
region («40 min). The real-time analysis capability was 
not available in 2010 but the data were also calibrated 
within AIPS following standard practice. 

As the experiment progressed we built up reference 
images, made by summing all previous epochs. These 
deeper images proved useful in the real-time search for 
transient sources (^5]). After the survey was completed, 
a final set of images was made separately for each yearly 
campaign using the data from the 11 epochs in 2008 and 
the five epochs in 2010. We also summed the 2008 and 
2010 deep images to create 16-epoch Master images for 
the entire experiment. In summary, there were three final 
image datasets; the Single epoch images, the Yearly im- 
ages (for 2008 and 2010 separately) and the final Master 
images made from all available data. 

Our final survey parameters are given in Table 2] The 
effective survey area was calculated using the full width 
at half-power (9.3' at 4.86 GHz; see however Appendix ICj) 
but the searches for transients in real time and for vari- 
ability were made over a larger area - out to the 15% 
response point of the primary beam (15' diameter). For 
the analysis, no correction was made for the primary 
beam attenuation, in order to maintain uniform noise 
statistics over the entire field. The synthesized beam 
and rms noise estimates for different epochs and differ- 
ent pointings varied by factors close to unity. The values 
in Table [3] are averages for each epoch over all point- 
ings, while Table |3] gives the mean rms values (over all 
fields) in the Master images and the 2008 and 2010 com- 
bined images. Throughout the paper we 
state explicitly if we use corrected or uncorrected fluxes. 
"Corrected fluxes" are corrected for beam attenuation 
and for the CLEAN bias by adding additional +0.3 mJy 
(e.g., Becker et al. 1995; Condon et al. 1998). In order to 
maintain uniform statistics we chose to search all images 
to the same depth, rather than compute a new threshold 
for each image individually. In practice, this led to some 
false positives for noisier than average epochs and fields 

10 http://www.aips.nrao.edu/ 



with bright point sources. Indeed, Table |3] indicates that 
the nosiest epochs contains larger number of sources. 

5. REAL-TIME TRANSIENT SEARCH 

We employed two distinct analysis strategies for tran- 
sient and variable source identification. The first, which 
is discussed in this section, was a real-time analysis. 
The main motivation was to rapidly identify any short- 
lived sources and mark them for immediate follow-up at 
other wavelengths. The second was a post-survey anal- 
ysis which was carried out after all the epochs had been 
observed. The main goals of this second phase were to 
carry out a more in-depth search for Bower et al. tran- 
sients (SjBJ) and to characterize the variability properties 
of the persistent source population ($7|). 

For the real-time identification the images ($4|) were 
searched visually for any new or strongly varying sources 
by comparing them with individual epochs, and by com- 
paring them with a reference image made by summing 
all previous epochs. Any candidate variable source which 
we identified was subject to a more detailed light curve 
and position fitting analysis before deciding to trigger 
radio, optical and/or X-ray follow-up observations. 

The followup visible light observations were carried out 
using the robotic Palomar 60" telescope (P60; Cenko et 
al. 2006) and the Keck-I 10-m telescope. The UV and 
X-ray observations were conducted by the Swift satellite 
(Gehrels et al. 2004). We note that prior to and during 
the VLA campaign we obtained visible light reference 
images for most of our fields using the P60 telescope. 

We identified two possible transients that were deemed 
interesting enough for multi-wavelength follow-up. How- 
ever, followup VLA observations and a careful post 
observing re-analysis (SJTJ) showed that these are not 
real transients. One source, J213438.01+414836.0, 
was a sidelobe artifact, while the second source, 
J230424.68+530414.7 is a long term variable that had 
crossed our single-epoch noise threshold on 2008 July 19 
and is clearly seen in the 2008 and 2010 deep coadds. 
We note that both sources were observed using the P60 
telescope about 2 hr and 1 hr after the radio observations 
were obtained, respectively. Furthermore Swift-XKF ob- 
servations of the first source were obtained about five 
days after it was found. These fast response observations 
demonstrate our near real time followup capabilities. 

6. POST SURVEY SOURCE CATALOG 

The next phase of our analysis occurred after the con- 
clusion of the observations. We generated source cata- 
logs in order to search for short-lived transients and to 
carry out a variability study of all identified sources. The 
AIPS task SAD (search and destroy) was used for source 
finding. 

We found that false positive sources came from one of 
two main reasons: slightly resolved sources and sidelobe 
contamination. Extended sources were identified by re- 
quiring that their integrated flux density was within a 
factor of two of their peak flux density. False sources 
created by scattered power from the snapshot sidelobe 
response was only a significant problem for the six fields 
with sources whose flux density exceeded 20mJy. We 
flagged any variables or transients from these fields for 
visual inspection. 
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TABLE 3 
Observing Epochs 



Epoch 


Date 
UTC 


Time Elapsed 
days 


(rms Noise) 

MJy 


Observed 
Fields 


Number 
of Sources 


Gain Corr. 


Cosmic error 

/o 


1 


2UUo Jul 15.4U 


0.00 


243 


111 


343 


1.029 


5.3 


2 


2008 Jul 18.39 


2.99 


181 


141 


155 


0.987 


1.6 


3 


nnno T.-i 1 n on 

2008 Jul 19.39 


3.99 


229 


141 


363 


1.033 


1.6 


■1 


nnno A.. in o o 

2008 Aug 10.33 


25.93 


178 


141 


166 


0.971 


1.0 


5 


2008 Aug 11.33 


26.93 


183 


141 


164 


0.980 


1.1 


6 


2008 Aug 14.32 


29.92 


171 


111 


162 


1.002 


0.4 


7 


2008 Aug 16.31 


31.91 


173 


139 


151 


1.023 


1.7 


8 


2008 Aug 18.29 


33.89 


178 


141 


155 


0.924 


1.7 


9 


2008 Aug 25.29 


40.89 


196 


111 


183 


1.016 


2.4 


10 


2008 Aug 28.28 


43.88 


178 


141 


157 


0.956 


3.5 


11 


2008 Aug 30.28 


45.88 


186 


141 


170 


1.033 


5.3 


12 


2010 Jul 16.42 


731.02 


105 


111 


216 


1.020 


1.9 


13 


2010 Jul 18.43 


733.03 


111 


134 


199 


1.022 


0.8 


14 


2010 Jul 22.35 


736.95 


108 


109 


158 


0.992 


0.6 


15 


2010 Jul 23.32 


737.92 


104 


111 


200 


1.008 


0.7 


16 


2010 Jul 25.31 


739.91 


116 


140 


217 


1.003 


2.1 


Note. 
which we 


— List of the 16 epochs. The dates indicate the observations mid-time. In practice, in all the instances in 
find that the cosmic error is smaller than 3% we replaced it by 3% (see ijol ■ 



TABLE 4 

Survey Parameters 



Property 


Value 


Frequency 


4.9 GHz 


Observing time 


48 hrs 


Survey Area 


2.66 deg 2 


Angular resolution 


15" 


Repeats 


16 


Timescales 


1 day-2 years 


No. of fields 


141 


mean exposure Time per field 


50 s 


mean rms per epoch (2008) 


190 /xjy 


mean rms per epoch (2010) 


109 A*Jy 


mean rms per 11 epochs (2008) 


72 /xjy 


mean rms per 5 epochs (2010) 


56 £tjy 


mean rms in Master images 


46.7/iJy 



Note. — 



Three catalogs were created. The first catalog is the 
"Single epoch catalog" , generated by running SAD on a 
15-arcmin diameter region for each single-epoch image 
individually (Table [5]). A second "Master catalog" (Ta- 
ble [5]) was created by running SAD on the final master 
images made from all available data (f6 epochs), while 
a third catalog, "Yearly catalog" , was generated on the 
yearly images (for 2008 and 2010 separately). 

The Master catalog has the merit of being able to iden- 
tify persistent radio sources approximately \/iV -times 
fainter than any individual epoch (where N = 16), but it 

is y/N -times less sensitive to a short-lived transient that 
might be identified in a single-epoch image. For our Sin- 
gle epoch catalogs we used a flux density cutoff of 1 mJy 
in 2008 and 0.76 mJy in 2010, and the number of sources 
in each epoch are given in Table [3l Our Master catalog 
consisted of 464 sources which are listed in Table [51 with 
a flux density cutoff of 0.28 mJy. The Yearly catalog had 
flux density cutoff of 0.5 mJy in 2008 and 0.35 mJy in 
2010. These cutoffs corresponds to about a 5-7er thresh- 
old, depending on the rms noise for individual fields. The 
electronic version of the Master catalog also contains the 
peak flux of each source. This was measured in the Sin- 



gle epoch images at the position of the sources found in 
the Master image. 

We used the Master catalog to perform a second order 
amplitude calibration that would tie together the flux 
density scale for all epochs. Normally self-calibration 
could be used to find additional gain variations within a 
radio observation but our survey was designed to avoid 
pointings with bright point sources. Our approach as- 
sumes that each VLA epoch (all the observations in each 
epoch were taken within 3 hours) shares the same "gain" 
correction, and we solved for these "nightly" gain cor- 
rections by fitting, using least squares minimization, the 
equation 

my = Zi + rhj, (1) 

where m»j is the "magnitude": — 2.5 log 10 fiji fij * s the 
peak specific flux of the j-th source in the i-th epoch, 
Zi is the gain correction for the i-th epoch (in units of 
magnitudes) and rhj is a nuisance parameter represent- 
ing the best fit mean magnitude of the j-th source. We 
note that, as explained in Appendix [A] magnitudes have 
convenient statistical properties. The final multiplicative 
gain corrections are 10 _0 ' 4z \ This method is described 
in detail in Appendix [AJ and the best fit multiplicative 
gain corrections are listed in Table 03 Similarly, we also 
derived the yearly gain corrections for the Yearly epochs. 
These gain corrections are 1.004 and 0.996 for 2008 and 
2010, respectively. 

The flux errors reported by SAD do not include any 
systematic error terms. Typically, VLA calibration is as- 
sumed to be good to a level of 3% or better (e.g., Condon 
et al. 1998). In order to check if some epochs are nosier 
we estimated the "cosmic errors", e cos , using the follow- 
ing scheme. We measured the standard deviation in the 
flux of the four phase calibrators observed on each night, 
after normalizing their flux by their mean flux over all 
the epochs taken at the same year. The cosmic errors 
estimated using this method are listed for each epoch in 
Table [3] This estimate is based on a small number of 
sources and these sources may be variable. Therefore, 
this should be regarded as a rough estimate. In some 
instances, the cosmic errors we estimated were smaller 
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TABLE 5 
Single epoch catalog 



Epoch 


Source 


Field Name 


J2000 RA 
deg 


J2000 Dec 
deg 


rcor 
J i/,p 

mjy 


runcor 

mjy 


CTp 

mjy 


f/fp 


Major 
ii 


Minor 
n 


PA 
deg 


Ac 
i 


1 


1 


1851 - 


1327 


282.955970 


-13.502335 


329.63 


25.39 


0.26 


1.03 


23.67 


13.67 


23.7 


2.89 


1 


2 


1851 - 


1327 


282.937502 


-13.432227 


3.29 


1.23 


0.26 


1.00 


23.49 


13.31 


29.7 


1.48 


1 


3 


1851 - 


1327 


282.855004 


-13.349517 


8.60 


1.12 


0.26 


0.64 


23.66 


8.48 


13.6 


8.30 


1 


1 


1853 - 


1233 


283.399554 


-12.561197 


7.59 


6.59 


0.26 


0.95 


21.88 


13.15 


20.3 


0.29 


1 


5 


1853 - 


1233 


283.474902 


-12.562245 


4.74 


3.10 


0.26 


0.89 


22.97 


11.74 


21.5 


4.13 



Note. — Catalog of 2953 sources detected in Single epochs. Columns description: Epoch is the epoch number (see Table [3t ■ 
and source is a serial source index in epoch. /^°p, /^p CDr , and a p are the corrected peak flux density, uncorrected peak flux density, 
and the error in the uncorrected peak flux density, respectively. Corrected flux are corrected for beam attenuation and the CLEAN 
bias, f/fp is the integrated flux divided by the peak flux. Major and Minor are the major and minor axes object size, while PA is 
the position angle of the major axis. Finally, Ac is the distance of the source from the beam center. This table is published in its 
entirety in the electronic edition of the Astrophysical Journal. A portion of the full table is shown here for guidance regarding its 
form and content. 



than 3% and in those cases we replaced the cosmic er- 
rors for these epochs by 3%. We note that if indeed the 
cosmic errors in some cases are smaller than 3%, then 
our strategy of adopting a larger cosmic errors may re- 
duce the number of variables found in our survey. For 
the Yearly catalogs, we used the mean cosmic error terms 
of the individual epochs in each year. These are 0.028 
and 0.012 for 2008 and 2010, respectively. 

Equipped with the gain corrections and an estimate 
for the cosmic errors we next corrected the flux measure- 
ments of all the sources using the gain correction factors 
and added in quadrature the cosmic errors to the peak 
flux errors. The new fluxes and errors were used in all 
the plots and the calculation of the light curves statistical 
properties. 

7. POST SURVEY TRANSIENT SEARCH 

Our final transient search utilized the catalogs pre- 
sented in fJH Specifically, we matched all the sources 
in the individual epochs to sources in the Master cat- 
alog using a 4" matching radius. Sources in individual 
epochs which do not have counterpart in the Master cat- 
alog are transient candidates. However, because we used 
a single flux density threshold, while the noise varies in 
each epoch and field, most of the faint sources are prob- 
ably noise artifacts. We used Figure [5] in order to choose 
a reasonable flux density limit for our transient search. 
This figure shows the number of detections, in the Single 
epoch catalogs, of sources which are detected in the Mas- 
ter catalog and their field was observed 16 times. This 
plot suggests that the probability that a faint source de- 
tected in the master image will be detected with uncor- 
rected flux density above 1.5 mjy, in only one epoch, is 
low. In fact it suggests that we could use an even lower 
threshold. However, given that this is based on small 
number statistics and given the variable quality of differ- 
ent images, we used an higher flux limit cutoff of 1.5 mjy. 

Following this analysis, we searched for sources de- 
tected in a single epoch that do not have a counter- 
part in the Master catalog, have uncorrected flux den- 
sity > 1.5 mjy, and a distance from beam center smaller 
than 4.65' (i.e., half power beam radius). In total we 
found 50 sources. However, a close inspection of these 
sources shows that most of them are not real. Of the 50 
candidates, 46 are in fields which contain sources brighter 
than lOmJy, and they are clearly the results of sidelobes. 
Of the remaining four candidates, three sources are also 
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Fig. 2. — Detection repeatability as a function of flux. The plot 
shows, for each source that its field was observed 16 times and that 
was detected in the Master catalog, the number of SAD detections 
in individual epochs. The number of detections is shown against 
the uncorrected peak flux of the source as measured in the Master 
image. Point sources are represented by filled circles, while open 
squares are for resolved sources. The vertical line shows the 1.5 mjy 
cut we used to define our completeness limit. 

most probably not real. One candidate was found next 
to a slightly resolved 5mJy source. A second object is 
a known 2 mjy source for which the centroid position 
puzzlingly shifted slightly in one of the epochs, a third 
candidate is a sidelobe seen in several epochs, and the 
fourth candidate is probably a real transient. This tran- 
sient candidate, J213622.04 + 415920.3 is described next. 

7.1. The transient candidate J213622.04 + 415920.3 

We found a single source, J213622.04 + 415920.3, that 
was detected only in the first epoch and may be a real 
transient. Given that this source was detected in the 
first epoch, before we constructed a reference image, it 
was not followed up in real time. The main properties of 
this transient candidate are summarized in Table [7] The 
peak flux measurements at the position of the source at 
all epochs show that the source was indeed visible only 
in the first epoch. Moreover, this source is not detected 
in the Master or Yearly images. The peak flux at this 
position in the Master image is 74 ± 42 /iJy, and in the 
2008 (2010) combined images is 161 ± 67 ^Jy (-14 ± 
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TABLE 6 

Master catalog 



Current search NVSS USNO 2MASS 

J2000 RA J2000 Dec /™ r /£ ncor a p N ohs N det x 2 StD/(/> x\ StD/(/) y Dist a a Dist Dist 
deg deg rnjy rnjy mjy " " " 



283.675661 


— 13.518385 


7.24 


2.48 


0.05 


15 


15 


50.23 


0.17 


2.54 





06 


1.0 


1.03 






283.682476 


— 13.444363 


2.52 


1.85 


0.05 


15 


15 


52.34 


0.18 


0.05 





02 


3.8 


0.17 






284.193460 


— 13.325061 


47.83 


32.01 


0.10 


14 


11 


70.30 


0.08 


0.48 





02 


0.9 


-0.07 




0.6 


284.387195 


— 12.213881 


11.03 


6.74 


0.06 


15 


15 


142.60 


0.14 


2.85 





05 


1.1 


-0.54 






283.976980 


— 12.165968 


93.19 


57.18 


0.13 


15 


15 


468.40 


0.18 


68.27 





25 


0.9 


0.08 






283.949467 


— 11.356239 


7.65 


4.69 


0.05 


15 


1 1 


53.20 


0.10 


0.84 





03 


2.9 


1.14 






no a r Aon r r* 

284.592955 


— 10.949209 


13.07 


10.66 


0.06 


15 


15 


48.94 


0.08 


3.98 





06 


4.2 


1.26 






284.497521 


— 10.459850 


2.70 


2.11 


0.05 


15 


14 


46.02 


0.14 


1.71 





06 


10.0 


0.83 






298.720898 


16.755257 


10.92 


8.39 


0.04 


16 


16 


118.76 


0.10 


5.02 





08 


1.0 


-0.80 






299.583716 


18.171732 


6.93 


5.82 


0.04 


16 


16 


143.31 


0.14 


15.03 





12 


2.3 


-0.20 






299.271523 


18.226378 


4.81 


3.65 


0.05 


16 


16 


75.27 


0.12 


0.67 





04 










299.353219 


18.473931 


41.09 


6.89 


0.05 


16 


16 


149.27 


0.13 


13.66 





13 


0.4 


0.14 






300.337088 


19.936532 


7.71 


5.50 


0.04 


16 


16 


55.62 


0.08 


0.12 





01 


2.2 


0.49 






300.030673 


20.034572 


23.81 


4.23 


0.04 


16 


16 


70.42 


0.11 


0.69 





02 


0.6 


0.39 










A 1A 


o.oo 


u.uo 


ID 


ID 


±0^.00 


die 




n 
u 


"1 7 
1 ( 










325.523312 


42. 023491 


12.85 


2.83 


0.04 


16 


16 


356.39 


0.20 


30.28 





20 


1.0 


0.88 






325.348791 


42.085546 


7.61 


2.23 


0.04 


16 


16 


148.57 


0.15 


8.95 





12 


0.3 


1.42 






326.519678 


Acs 1 a i ore? 

42.141397 


30.49 


8.09 


0.04 


16 


16 


54.95 


0.06 


2.98 





05 


0.5 


1.09 






325.966278 


42.8844o3 


7.55 


2.34 


0.05 


16 


16 


46.73 


0.15 


0.45 





02 










326.391328 


43.460876 


7.24 


4.02 


0.04 


16 


16 


46.97 


0.10 


0.32 





01 


1.8 


1.00 






341.141889 


51.541980 


15.60 


12.43 


0.05 


16 


16 


160.78 


0.10 


0.53 





03 


1.6 


0.54 


0.3 




342.654997 


52.100799 


24.51 


13.77 


0.05 


16 


16 


195.40 


0.11 


18.74 





13 


1.8 


0.28 


0.6 


0.5 


341.820302 


52.137525 


4.11 


2.64 


0.04 


16 


16 


49.55 


0.12 


0.43 





03 


1.9 


0.66 


0.4 


1.1 


343.612547 


co oonm n 

52.339910 


9.04 


7.77 


0.05 


16 
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Note. — Catalog of the 464 sources detected in the Master images and their properties. We note that SAD detected additional sources which arc 
not listed here. These sources were identified as noise artifacts by subsequent inspection of the images and were removed from the catalog. Columns 
description: /£°p, f^ U p OT , and cr p arc described in Table[5] V b s is the number of epochs in which the source position was observed, while N^ ct is the 
number of detections in the Single epoch images. Subscript L 'Y" in x 2 and StD/(/) indicates that these values are calculated for the Yearly epochs. 
Dist is the distance between the radio position and the NVSS, USNO-B1 and 2MASS (Skrutskic ct al. 2006) nearest counterparts. A portion of the 
full table, containing the 30 Single epoch variable sources ( EI8.H and the ten Yearly variable sources ( arc shown here for guidance regarding 
its form and content. The two variable lists arc separated by horizontal line. This table is published in its entirety in the electronic edition of the 
Astrophysical Journal. The electronic tabic also contains additional columns as the Field name, sidclobes flag, integrated flux divided by peak flux, 
major axis, minor axis, and position angle of the sources, distance from the beam center, Vr. and Vp for the yearly fluxes, the NVSS flux, the B and i? 
magnitudes of the USNO-B counterparts, J, H and ./^-magnitudes of the 2MASS counterparts and all the 16 peak fluxes and errors measurements in 
the Single epochs and the peak flux measurements in the Yearly images. We note that the distance threshold for USNO-B and 2MASS counterparts 
was set to 2", and 15" for NVSS. 

a The spectral power-law slope, defined by /j, cc u a , as measured from the NVSS 1.4 GHz specific flux and our 4.9 GHz corrected flux density. 



49 /i Jy) . In order to test if the source is variable 

during the 46 s VLA integration, we split the data into 
two 23-s images. We found that the flux of the source 
did not change significantly between the first and second 
parts of the exposure. 

The field of J213622.04 + 415920.3 was observed with 
the P60 telescope in i-band about two days after the 
transient was detected. This image is shown in Figure [31 
We find two sources near the transient position. One 
source is found = 3.2" North-East of the transient loca- 
tion and has an i-band magnitude of 20.98 ± 0.23. The 
other is found = 2.3" South of the transient position 
and has a magnitude of 19.46 ± 0.20. The magnitudes 



are calibrated relative to USNO-B1.0 I-band magnitude 
(Monet et al. 2003). Given the relatively large angular 
distances between the transient position and the nearest 
visible light sources, they are probably not associated. 
We note that given the stellar density in this region the 
probability to find a source within a 2.4" radius from a 
random position is about 15%. Finally, we do not find 
any counterpart to this transient in the SIMBAD, NED 
or HEASARC databases. 

The search method described in the beginning of this 
section may miss transients which are bright enough to 
be present in the Master catalog. Therefore, we also 
searched for sources which are detected in the Master 
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TABLE 7 

Properties of the transient candidate 
J213622.04 + 415920.3 



Property 



Value 



Field name 
R.A. (J2000.0) 
Dec. (J2000.0) 
Detection date 
Uncorrected peak flux 
Corrected peak flux 



2136+4158 
21 h 36 m 22. s 04± 1.2" 
+41°59'20."3± 1.2" 
2008 Jul 15.4147 
1.61 ± 0.28 mjy 
2.36 ± 0.41 mjy 



Distance from beam center 2.77' 

Note. — We note that the coordinates in this ta- 
ble are based on Gaussian fit. The coordinates of 
this source in the Single epoch Table (Table O wore 
derived using SAD and they are somewhat different 
RA= 21 h 36 m 21. s 965, Dcc= +41°59 / 21. // 52 (J2000.0). 
Note that the flux is corrected also for the CLEAN bias. 



J21 3622.04+41 5920.3 
17-Jul-2008 P60/i-band 



• 



Q 




Fig. 3.— P60 i-band image of the field of J213622. 04 + 415920.3. 
The image, with exposure time of 540 s, was taken on 2008 Jul 
17.45. A 2.4" radius circle marks the transient candidate location. 
This radius roughly corresponds to the 2-cr error circle. 

catalog and detected in only one of the Single epoch cat- 
alogs with uncorrected flux density above 1.5 mjy. No 
such sources were found. 

8. VARIABILITY ANALYSIS 

We investigated the variability of all the sources in the 
Master catalog using their peak flux densities measured 
in the Single epoch images f §8.1|) . and the Yearly images 
( 38.2[) . We used various statistics to assess the variability 
and its significance, which we list in Table [51 For each 
source we calculated its standard deviation (StD) over 
all the flux measurements, the StD over the mean flux 
(StD/(/)) and the x 2 given by 



N 

E 



(fi-(f)) 2 

°f + (/i£cos,i) 2 ' 



(2) 



where N is the number of measurements which is 2 for 
the Yearly catalogs and 16 for the Single epoch catalogs, 
i is the epoch index, is the gain corrected peak flux 
in the i-th epoch, <7j is its associated error, e cos is the 
cosmic error, and (/) is the mean flux of the source over 



all the epochs. Note that in Table HI x 2 is measured over 
the individual epochs, while \\ is measured over the two 
yearly epochs. In the first case, the number of degrees of 
freedom (dof) is 15, while in the second case it is one. 

Some previous surveys have defined "strong variables" 
as exceeding some pre-defined variability measure. There 
are many definitions of fractional variability in the litera- 
ture and for comparison with other surveys we list in the 
electronic version of Table [6] the two following indicators 



max{/. t } 
min{/j} : 



and 



max{/,} + min{/J ' 



(3) 



(4) 



We note that these two quantities are related through 
Vf = (Vr — 1)/(Vr + 1). However, as neither of these 
indicators account for measurement errors, they cannot 
be used reliably on their own at low fluxes where the 
measurement errors are very large. Most importantly, 
estimators involving the min or max functions strongly 
depend on the number of measurements. This fact com- 
plicates direct comparison between different surveys. For 
example, a light curve whose fluxes are drawn from a log 
normal random distribution, with N cp = 3 survey, and 
StD/(/) = 0.5 has (Vr) = 2.44, while for N cp = 16 it 
will have (Vr) = 5.62. Therefore, although Becker et 
al. (2010) and Taylor and Gregory (1983) defined strong 
variables identically, i.e. Vr >3 (Vf > 0.5), any di- 
rect comparison of these two surveys is difficult since the 
number of epochs in these surveys were 3 and 16, respec- 
tively. In the future, in order to cope with this problem 
we suggest to use the StD / (/) as a fractional variability 
estimator. Moreover in Appendix [B] we provide a con- 
version table for Vr as a function of AT ep and StD/(/). 

8.1. Short time scale variability 

In order to explore radio variability on short time scales 
(e.g., days to weeks), we constructed 16 epochs light 
curves for all sources with peak flux densities larger than 
1.5 mjy (w 6er). Figure 0] shows the StD/(/) vs. the x 2 
of all sources in the Master catalog. This figure suggests 
that a large fraction of at least the bright radio sources 
with flux densities larger than about lOmJy are variables 
at the level of £ 5%. 

In total, we find that 30% (30 out of 98) of the sources 
in our survey, which are brighter than 1.5 mjy, are vari- 
able (at the 4-cr leveO). The light curves of these 30 
variable sources are presented in Figure [5J and their flux 
measurements and basic properties are listed in Table [6j 
This is considerably larger than the fraction of variables 
reported in some of the other "blind" surveys listed in 
Table[H(e.g. Gregory & Taylor 1986; de Vries et al. 2004; 
Becker et al. 2010). A possible explanation for this ap- 
parent discrepancy is that the sources in these surveys 
were extracted from mosaic images in which each point 
in the survey footprint was observed multiple times dur- 
ing several days. Therefore, the fluxes they reported are 
average fluxes over several days time scales (column St 

11 Assuming Gaussian noise, 4o" corresponds to a probability of 
= 1/15, 000 while the number of measurements in our experiment 
(number of epochs multiplied by the number of sources) is 7400. 
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Fig. 4. — StD/(/} vs. the x 2 , where the uncorrected peak flux 
of the sources (in the Master catalog) is marked by symbol size. 
The dashed line corresponds to x 2 > 45.5, which corresponds to 
4cr assuming 15 degrees of freedom. 

in Table [1]) . Such measurements will tend to average out 
variability on time scales which are shorter than St. In 
order to test this hypothesis we carried out a structure 
function analysis. 

We calculated the mean discrete auto-correlation func- 
tion, C(t), of all the 30 variable sources, as a function 
of the time lag r. We first normalized each source light 
curve by subtracting its mean and than dividing it by 
its (original) mean. We treated all these light curves 
as a single light curve by concatenating them with gaps, 
which are larger than the time span of each light curve, in 
between light curves. Then we followed the prescription 
of Edelson & Krolik (1988) for calculating the discrete 
auto-correlation function. The errors were calculated us- 
ing a bootstrap technique with 100 realizations for the 
measurements in each time lag (e.g., Efron 1982; Efron 
& Tibshirani 1993). The mean auto-correlation func- 
tion is presented in Figure [6] (black circles) . The auto- 
correlation at lag "zero" is not one. This is because at lag 
zero we used a lag window of to +2 days. Therefore it 
does not contain only zero lag data. The auto-correlation 
function reaches zero correlation at r « 10 days. 

Next, we calculated the structure function, SF(t), of 
these light curves defined by 

2S' 2 Ar'M 

SFM = \/2ff 2 (l - CM) ± - U= , (5) 

y ' V V L U 2 V /2S 2 (1-C[r]) 

where S is the standard deviation of the normalized light 
curves and AC[r] is the bootstrap error in the discrete 
auto-correlation function. The second term on the right 
hand side of Eq. [5] represents the error in the structure 
function. The structure function is presented in Figure[7J 
Also shown in Figures [6] and [7] are the auto-correlation 
and structure functions, respectively, for "non-variable" 
sources (gray symbols). In this context our selection cri- 
teria for non- variable sources are specific flux larger than 
2mJy and \ 2 < 25.3. This \ 2 value corresponds to 2-a 
confidence, assuming 15 degrees of freedom. 

The structure function of the variable sources, after 
subtracting the non- variable source structure function is 



rising rapidly from zero to sa 0.05 on a time scale of the 
order of one day and than rises to a level of ~ 0.12 at lags 
of m 10 days at which it stays roughly constant. Though, 
we cannot rule out it is slowly rising on r > 10 days time 
scales. 

This analysis suggests that a large component of the 
variability happens on time scales shorter than about 
one day. A plausible explanation is that the short time 
scale variability is due to refractive scintillations and it 
is discussed in gS21 (e.g., Rickett 1990). This level of 
variability was likely missed by some previous surveys 
(Table [T]) due either to their choice of frequency (v), ca- 
dence (N ep ) or the observing time span (St). On the 
other hand, 5 GHz surveys of flat spectrum active galac- 
tic nuclei (AGN) find that their majority shows signif- 
icant variability, and that the number of these variable 
sources increase at low Galactic latitudes (e.g., Spangler 
et al. 1989; Ghosh & Rao 1992; Gaensler & Huntstead 
2000; Ofek & Frail 2011). The source population in the 
flux density range of our survey is known to be domi- 
nated by AGN and a significant fraction (~ 50%) of these 
are compact, flat-spectrum AGN and hence expected to 
show short-term flux density variations (de Zotti et al. 
2010). 

8.2. Variability on years time scale 

By comparing our 2008 and 2010 catalogs we were 
able to probe the variability of sources brighter than 
J> 0.5 mJy on two year time scale. However, this com- 
parison is limited to only two epochs. Each of the two 
epochs in our two year time-scale variability analysis is 
composed of multiple epochs taken (« St) days to weeks 
apart. Therefore, in this analysis short-term variability 
is averaged out. For example, in 5 GHz, variability due 
to scintillations will typically have time scale which is 
shorter than a few days. 

Figured] shows a comparison of the 2008 and 2010 peak 
flux densities for all radio sources in the Master catalog. 
The dashed lines represent the mean noise 4 and 8-a con- 
fidence variability contours. Based on the x\ (Table O 
we find ten sources (2% of all the sources) which have 
variability with confidence level larger than about 4a 
(i.e., Xy > 16; assuming one degree of freedom). These 
variable candidates are listed in Table [6] (below the hor- 
izontal line). We note that the x\ calculation takes into 
account the gain correction factors and the cosmic errors 
described in Sj6] We visually inspected the images of all 
ten variable candidates and verified they are not sidelobe 
artifacts. 

J213626.36 + 413926.6 is the only significant strong 
variable (i.e. with x 2 /dof > 16/1 and V-p > 0.5) in 
our two-year comparison. In 2010 it has a peak flux 
density (uncorrected for primary beam attenuation) of 
627 ± 50/iJy (12.5er), while in 2008 there is a nominal 
detection of the source at this position with a peak flux 
density of 197 ± 75/xJy (2.6a). The source flux was well 
below the detection threshold for every one of the 11 
epochs in 2008 but it is visible in all five epochs in 2010. 
We note that this source is one of the two faint sources 
seen above the dashed line on the right side of Figure |H 
For this reason we do not classify this as a transient 
source but it appears instead to be a persistent radio 
source that has tripled its flux density over a two year 
interval. There is no known cataloged source at this po- 
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Fig. 5. — Light curves of the 30 variables which have \ 2 > 45.5. For scaling purposes we show only the 11 observations taken during 
2008. In each panel we give the J2000.0 Right Ascension and Declination of the source. The individual flux measurements are given in the 
electronic version of Table [6] 



sition in the NVSS catalog (Condon et al. 1998) nor in 
the SIMBAD or NED databases. 

Our Master catalog contains 317 sources brighter than 
J213626.36 + 413926.6. Therefore, we roughly estimate 
that the fraction of strong variables with two epochs 
V F > 0.5 is 0.32t°;£?%. However, we note that this mea- 
surement is based on averaging out variability on time 
scales shorter than a few weeks. 



9. DISCUSSION 

We present a 16 epoch, Galactic plane survey, for radio 
transients and variables at 5 GHz. We detected one pos- 
sible transient and many variable sources. The transient 
areal density and rate based on this single detection are 
derived in £19.11 In §9.21 we discuss our variability study 
and compare it with previous surveys. 
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Fig. 6. — The mean discrete auto-correlation function of all the 
30 variable sources (black circles). The gray circles show the same 
but for all the non- variable sources (see text). Each light curve was 
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Fig. 7. — The mean structure function of the normalized light 
curves as a function of lag, r. Symbols like in Figure [6] Fol- 
lowing Lovell et al. (2008) we fit the structure function (after 
subtracting the non-variable sources structure function) with the 
function m s T/(r + T c i lar ). The heavy gray line represent the best 
fit function, and the best fit parameters are m B = 0.119 ± 0.003, 
Tchar = 0.83±g:lJ day ( X 2 /dof = 26.6/6). 



9.1. Transient areal density and rate 

The analysis of our data revealed a single radio tran- 
sient candidate. The area encompassed within a single 
half power beam (radius of 7'hp = 4.65') in which we 
searched for transients is 0.00601 deg 2 and the area we 
targeted within the half power radius of all 141 fields is 
2.66 deg 2 . Given that each field was observed on aver- 
age 15.70 times (sum of column four in Table [2] divided 
by the number of fields), the total area covered by our 
survey over all the epochs is 41.2 deg 2 . Our survey used 
an uncorrected flux limit of 1.5 mJy for transients. How- 
ever, the sensitivity within the field of view of a single 



beam imaging is not uniform, and degrades by a factor 
of two at the half power radius. In order to calculate the 
transient areal density from these parameters we need to 
assume something about the source number count func- 
tion. We parametrize the source number count function 
as a power law of the form 



*(> /) = «o(///o)- 



(6) 



where / is specific flux, n(> f) is the sky surface density 
of sources brighter than /, no is the sky surface density 
of sources brighter than fo, and a is the power law index 
of the source number count function. It is well known 
that for homogeneous source distribution in an Euclidean 
universe and arbitrary luminosity function a = 3/2. In 
Appendix [C] we derive a simple relation for the number 
of sources that are expected to be detected in a beam 
with a power sensitivity that falls like a Gaussian as a 
function of fo, kq, a, 7'hp, the search radius r max , and 
the specific flux limit at the beam center / m i n .o- 

Based on Equation lC51 and assuming a — 3/2, we find 
that the transient areal density at 1.8 mJy (corrected for 
the CLEAN bias) is 



k(> 1.8 mJy) 



non+0.13,+0.18 , -2 
u - uo3 -0. 032,-0. 038 ue & ' 



(7) 



where the errors correspond to one and two sigma con- 
fidence intervals, calculated using the prescription of 
Gehrels (1986). If our detected transient is not real then 
our survey poses a 95% confidence upper limit on the 
transient rate of k(> 1.8 mJy) < 0.15 deg -2 . 

Translation of our areal density to transient rate de- 
pends on the transient duration i^ur and it is 



5ft(> 1.8 mJy) = (28; 



-65, + 132 
-23,-27 



t,h. 



0.5 day/ 



I deg 2 yr 1 . 



( 8 ) 

Note that this translation is correct only if tdur is smaller 
than the time between epochs. 

Figure Q] presents a summary of the radio transient 
and persistent source areal density as observed by vari- 
ous searches and at different frequencies. This figure is 
largely based on Table Q] and the areal density reported 
here. As shown in this figure, the areal density derived in 
this work is consistent with the expectation based on the 
Bower et al. (2007) transient sky surface density. More- 
over, it is roughly consistent with the sky surface density 
of the Nasu survey transients. We note that the Nasu 
sky surface density is based on our limited knowledge of 
this project (see $2]). This comparison assumes that the 
transient areal density on the celestial sphere is uniform. 

This figure implies that the areal density of radio tran- 
sients in the sky is roughly 2-3 orders of magnitude below 
the persistent radio source sky surface density. This is in 
contrast to the visible light sky in which the fraction of 
transients (excluding solar system minor planets) among 
persistent sources is roughly 10~ 4 down to the limiting 
magnitude of surveys like the Palomar Transient Factory 
(Law et al. 2010; Rau et al. 2010). 

It is interesting to compare this figure with some re- 
cent predictions. Nakar & Piran (2011) predict that com- 
pact binary mergers, regardless whether they are associ- 
ated with short-duration gamma-ray bursts (e.g., Nakar 
2007), will produce radio afterglows with a duration of 
several months. They suggest that the two-months long 
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Fig. 8. — A comparison of the 2008 and 2010 peak flux densities for all radio sources in the Master catalog. Equal fluxes are indicated 
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a flux dependent component which is given by the cosmic errors discussed in ij6] The average rms noise for the summed 2008 images was 
72/iJybeam — 1 and for 2010 images it was 56 fijj beam - 1 . The shape and color coding represent the x 2 (one dof) of individual objects 
based on the actual noise rather than the average noise (see legend). Open triangles mark limits were the direction of the limit is given by 
the direction of the tip of the triangle (bottom or left). 



radio transient RT19870422 detected by Bower et al. 
(2007) may be a binary merger radio afterglow. More- 
over, they find that the rate inferred from this event 
is consistent with the predicted rate of binary merger 
events. 

Giannios & Metzger (2011) suggested that tidal flare 
events may produce radio transients with durations of 
month to years, with a 5 GHz peak flux of 1 mJy at a 
distance of 1 Gpc (i.e., z = 0.2). The total comoving 
volumePl enclosed within a luminosity distance of 1 Gpc 
is 2.4 x 10 9 Mpc 3 or 5.9 x 10 4 Mpc 3 deg~ 2 . Bower et 
al. (2007) did not find any transients with a duration 
of two months which are associated with the nucleus of 
a galaxy. This is translated to a 95% confidence upper 
limit on the rate of radio tidal flare events, brighter than 
lmJy, of ~ 0.1 deg -2 yr _1 . Therefore, in the context 
of the Giannios & Metzger (2011) predictions, we can 
put an upper limit of ~ 7 x 10~ 6 Mpc 3 yr" 1 on the rate 
of tidal flare event radio afterglows. This is in rough 
agreement with the predicted tidal flares rate (Magorrian 

12 Assuming WMAP fifth year cosmological parameters (Ko- 
matsu et al. 2009). 



& Tremaine 1999; Wang & Merritt 2004). 

Finally, we note that our detection rate is consistent 
with an old NS origin as suggested in Ofek et al. (2010), 
and with their expected surface density at low Galactic 
latitude based on the Ofek (2009) simulations. 

9.2. Comparison of variability with previous surveys 

We analyzed the variability of the sources detected in 
our survey on days-weeks time scales and two years pe- 
riod. On short time scales we find that considerable frac- 
tion of the bright point sources are variable. At the 10- 
100 mJy range it seems that more than half the sources 
are variable on some level (> 4a). Furthermore, we find 
that 30% of the sources in our survey, which are brighter 
than w 1.5 mJy, are variable (at the 4-a level). This 
is considerably higher than the fraction of variables re- 
ported in some other surveys (e.g. Gregory & Taylor 
1986; de Vries et al. 2004; Becker et al. 2010). We sug- 
gest that a possible explanation for this apparent discrep- 
ancy is that the sources in these surveys were extracted 
in a way that washed out short time scale variability (see 
38.ip . This is supported by the fact that our two years 
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time scale variability study, in which each epoch is com- 
posed by averaging multiple observations, shows smaller 
fraction of variables. Moreover, our structure function 
analysis shows that a big fraction of the variability com- 
ponent happens on time scales shorter than about a few 
days. We note that these fast variation of radio sources is 
known for a long time, and was found by Heeschen (1982; 
1984). Moreover, a large fraction of variable sources was 
previously reported by some other efforts, which did not 
average out short time scales variability (e.g., Lovell et 
al. 2008). 

We speculate that the fast rise of the structure function 
on ss 10 day time scales is due to scintillations in the in- 
terstellar medium (ISM). These time scales are consistent 
with those expected theoretically from refractive scintil- 
lations in 5 GHz (e.g., Blandford et al. 1986; Hjellming 
& Narayan 1986). Moreover, similar rise times were re- 
ported by other efforts (e.g., Qian et al. 1995). Unlike 
diffractive scintillation which may produce strong vari- 
ability (StD/(/) of 80%; e.g., Goodman 1997; Frail, Wax- 
man & Kulkarni 2000), refractive scintillations can easily 
explain the observed amplitude of sa 13%. We note that 
a comparison of models with observations suggests that 
most of the radio source variability below 5 GHz is due 
to scintillations, while above 5 GHz there is an intrinsic 
variability component (e.g., Hughes, Aller & Aller 1992; 
Mitchell et al. 1994; Qian et al. 1995) Therefore, we can 
not rule out the possibility that some of the variability 
we detected in our survey is intrinsic to the sources. 

After averaging out variations on days time scale, our 
two year variability analysis indicates that about 0.3% of 
the sources above 0.5 mJy are strong variables (Vf > 0.5 
for 7V ep = 2), and that only a small fraction s» 3% of the 
sources are variables at some level. This finding supports 
the hypothesis that the main reason for low-amplitude 
radio variability at 5 GHz is due to scintillations. 

It is well known that the fraction of radio variable 
sources increases toward the Galactic plane (e.g., Span- 
gler et al. 1989; Ghosh & Rao 1992; Gaensler & Hun- 
stead 2000; Lovell et al. 2008; Becker et al. 2010; Ofek 
& Frail 2011), plausibly due to Galactic scintillations. 
However, there are some claims as yet unconfirmed that 
the number of intrinsically strong variables, at 5 GHz, 
increases toward the Galactic plane (Becker et al. 2010). 
Specifically, Becker et al. (2010) suggested that there is 
a separate Galactic population of strong variables. As 
noted in Sj2j they found more than half of their variables, 
on timescales of years, varied by more than 50% in the 1- 
100 mJy flux density range. Moreover, they found that 
these strong variables were concentrated at low Galac- 
tic latitudes and toward the inner Galaxy. In contrast, 
we find a much smaller fraction of strong variables (see 
§8.2j) . However, Becker et al. (2010) observed sources 
within one degree of the Galactic plane, while our sur- 



vey sampled Galactic latitudes |6| = 6-8°. 

10. SUMMARY 

We present a VLA 5 GHz survey to search for radio 
transients and explore radio variability in the Galactic 
plane. Our survey represent the first attempt to dis- 
cover radio transients in near real time and initiate multi- 
wavelength followup. Our real time search identified two 
possible transients. However, followup observations and 
our post-survey analysis showed that these candidates 
are not transient sources. Nevertheless, in one case, we 
were able to initiate visible light observations of the tran- 
sient candidate field only one hour after the candidate 
was detected. 

Our post survey analysis reveals one possible transient 
source detected at the 5.8-cr level. Our P60 images of 
this transient field, taken two days after the transient 
detection, do not reveal any visible light source brighter 
than i-band magnitude of 21 associated with the tran- 
sient within 2". The transient has a time scale longer 
than 1 min. However, we cannot put an upper limit on 
its duration since it was detected on the first epoch of our 
survey. Based on this single detection we find a transients 
brighter than 1.8 mJy areal density of f=s 0.04 deg -2 . The 
transient surface density found in this paper is compared 
with other surveys in Figure [T] Our transient areal den- 
sity is consistent with the one reported by Bower et al. 
(2007), corrected for the flux limit. This is also roughly 
consistent, up to a possible spectral correction factor, 
with the rates reported by Levinson et al. (2002) and 
Kida et al. (2007). Finally, based on existing evidence, 
we cannot rule out the hypothesis that these transients 
are originating from Galactic isolated old NSs (Ofek et 
al. 2010). 

Finally, we present a comprehensive variability analy- 
sis of our data, with emphasis on proper calibration of 
the data and estimating systematic noise. Our findings 
suggest that short time scale variability among 5 GHz 
point sources is common. In fact above 1.5 mJy at least 
30% of the point sources are variable with variability ex- 
ceeding our 4-er detection level. This is consistent with 
the Lovell et al. (2008) results, and is plausibly explained 
by refractive scintillations in the ISM. 
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APPENDIX 

A. ESTIMATE OF THE GAIN CORRECTION FACTORS 

In order to construct the best light curves of the sources we need to remove any systematic factors influencing the 
measurements. A way to do this is to use the fact that the actual light curves of many sources are not correlated. 
Therefore, in the absence of systematic factors affecting the measurements, the scatter in the average light curve should 
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be minimized. In order to minimize the scatter in the average light curve of our sources we multiplied the fluxes at each 
epoch of observation (taken during ~ 3 hours) by a "gain" correction factor, such that the residuals in all light curves, 
compared with a constant light curves, will be minimized. This problem, is similar to producing relative photometry 
light curves in optical astronomy. 

We used a linear least square minimization technique. Using this method we are solving for the best zero point 
normalization (per epoch), and the best "mean" flux of each source, that minimize the global scatter in all the light 
curves. This technique was already introduced by Honeycutt (1992), but here we write it in a more easy to use form 
and we also add linear set of constraints for simultaneous absolute calibration. 

We work in a "magnitude" system (i.e., m = — 2.51og 10 /). The advantage of the magnitude system is that the 
logarithm of a quantity has the property of making error distribution more symmetrical, and this may be somewhat 
important for faint sources. The basic idea of this technique is to simultaneously solve the following set of equations 
in the least square sense: 

niij = Zi + fhj, (Al) 

where m^ is a p x q matrix that contains all the measured ("instrumental") magnitudes, i is the epoch index (p 
epochs), and j is the source index (q sources). Here, fhj is the mean magnitude of the j-th source and is the zero 
point of the i-th epoch. We note that fhj and Zi are free parameters. In case we have error measurements, let <r£? 
be the respective errors in the instrumental magnitudes. In some cases, we may have additional constraints like the 
calibrated magnitudes, Mj, (and respective errors, <rj^) of some or all the sources. This additional information can be 
used as constraints on the system of linear equations. Using these constraints our output magnitudes will be calibrated 
in respect to a set of reference sources. 

Given m,j and <r™, and the optional Mj and crj^ we would like to find the properly weighted best fit free parameters 
Zi and fhj. Let m be a vector of the observable quantities (reorganization of m^), and a the respective vector of 
errors in these quantities obtained by re-arranging the matrices of instrumental and calibrated magnitudes (i.e., my 
and aVj): 

mu 
mi2 



'la- 




in. = 



M n 



a = 



'21 



'2q 



'pi 
T p2 



(A2) 



Note that the elements below the double horizontal line are optional elements, needed only if we want to simultaneously 
apply a magnitude calibration. 
Next, we can define a vector of free parameters we would like to fit: 



P=[zi z 2 ... 



mi m-2 



m a 



(A3) 



where the superscript T indicate a transpose operator. In that case the design matrix H, which satisfy m = HP is 
easy to construct: 



H 



■ T qxp 
.7=1 


jqxq ' 


j-qxp 


jqxq 




jqxq 


jqxp 


jqxq 


.7=9 


Qqxp 


jqxq 



(A4) 
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where 



is a q x p matrix in which the A:-th column contains ones, while the rest of the elements are zeros, is 



a q x q identity matrix, and qxp is a q x p zeros matrix. 

Note that (again) the lower block of the matrix H, separated by two horizontal lines, is an optional section that 
is used for the magnitude zero-point calibration. This additional section acts like constraints on the system of linear 
equations. 

The rank of the matrix H without the lower block is p + q — 1. This is because without the calibration block there 
is an arbitrariness in adding a zero-point to each epoch. Adding the calibration block (or part of it) fixes this problem 
and in that case the rank of H is p + q. In cases when a given source does not appear in a specific epoch, we simply 
have to remove the appropriate row in H, rh and a . 

In order to find the best fit parameters P, and their respective errors a p , we need to find P that minimizes the x 2 : 



described in many textbooks (e.g., Press et al. 1992; for a tutorial see Gould 2003). 

The design matrix, H, even without the magnitude calibration part, is an {p x q) x (p + q) matrix. For many 
problems, the matrix H, may be huge and requires a lot of computer memory. However, H is highly sparse, and 
therefore sparse matrix utilities can be used if needed. Alternatively, this can be solved using the conjugate gradient 
methocF^I. 

We note that in practice this method may be applied iteratively. After the first iteration, we can check the Xj f° r 

each source and the % 2 for each epoch. Then, we can remove sources with large value of x 2 , and/or we can add a 

"cosmic" errors to all the measurements in an epoch with large x? • After we construct the new H and Y we may 
apply the inversion again. 

We note that typically, in addition to a™ (the errors associated with the individual sources) there are additional 
errors (e.g., calibration errors in radio astronomy and flat fielding errors in optical astronomy). Ignoring these errors 
is not recommended since it will give over weight to sources with small errors. A solution to this problem is, again, to 
apply this method in iterations. After the first iteration it is possible to estimate the cosmic error term (based on the 
residuals from the best fit), and to add these cosmic errors to the instrumental errors in the second iteration. 

Finally, using this method additional de-trending is possible. For example, one can add additional columns to the 
design matrix (and corresponding additional terms to the vector of free parameters) that represent changes in the zero 
point as a function of additional parameters. For example, in radio astronomy, this may be a variation in the zero 
point as a function of the distance from the beam center, and in optical astronomy this may be an airmass-color term, 
positional terms, color terms affecting different instruments and more. 



The Vr and Vp defined in Equations [3HH are sensitive to the number of measurements. This complicates the 
comparison between surveys with different number of epochs. In order to demonstrate this and to provide with a mean 
for roughly convert these variability indicators between different surveys, here we calculate the expectation value for 
Vr as a function of the number of epochs in a survey, -/V ep , and the StD/ (/) of a source light curve. We note that Vr 
and Vf are exchangeable. 

In order to calculate this conversion we performed the following simulations. We generated random light curves with 
iV e p points and which are drawn from a log- normal standard deviation, SiogN- For each value of iV e p and SiogN 10 
light curves were generated and StD/(/) and (Vr) were calculated. In Table |8] we list the Vr expectation values as a 
function of N ep (rows) and StD/(/) (columns). Below the StD/(/) we give also the appropriate Si og N- 



Typically, the sensitivity of a radio telescope is not uniform across its field of view, and depends on the radial angular 
distance from the beam center. In order to convert areal density, kq, of sources brighter than flux density /o, to the 
expected number of detectable events by a radio telescope we need to take into account the beam pattern and the 
source number count function. 

We parametrize the cumulative density of events as a function of flux as a power law 



where n(f)df is the number of sources per flux density interval, and a is the power-law index of the source number 
count function. For a uniform density population with arbitrary luminosity function in an Euclidean universe a — 3/2. 

The number of sources that can be detected in a single beam in a single epoch up to an angular distance r max from 
the beam center is 




B. Vr AND Vf STATISTICS 



C. ESTIMATE OF THE AREAL DENSITY IN A NON UNIFORM BEAM 




(CI) 
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TABLE 8 
Translation of (Vr) to StD/(/> 
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'-'logN 
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Note. 



27rrre [/minM//o] 



(C2) 



where r is the distance from the beam center, / is the flux density, and / m in(?") is the detection threshold as a function 
of angular distance r. For convenient we will assume that the beam pattern is Gaussian so that 



f • (r) - f • ne + r2ln2 /( r Hp) 

where rnp is the half width at half poweJ^I. and / m in,o is the detection limit at the beam center (i.e., r = 0). 
In the case of a = 3/2 and a Gaussian beam pattern, the integral in Equation IC2I has an analytic solution 



N h (a = 3/2) 



27r« r| p 
31n2 



min,0 c +r 2 ln2/r^ p 



fo 



-3/2 



(C3) 



(C4) 



For other values of a this integral can be evaluated numerically. Re-arranging Equation IC4| the surface density, 
assuming a — 3/2, is given by 



3iV 6 ln2// mm , x3/2 



2nr 2 HP 



fo 



_ e -3<„ ax ln2/(2r^ p ) 



(C5) 



Finally, the ratio between 7rr^ ax and N(a = 3/2;r max ) gives the correction factor of 1.61 we used in 
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